package com.mg.plugin.agrometeorology.common;

public class PenmamCalculator {

    /**
     * ta平均温度，wd风速，day日序号，y总日数，hd相对湿度，rizhao总辐射，H高度，LAT纬度
     *
     * @param ta
     * @param tmax
     * @param tmin
     * @param wd
     * @param day
     * @param y 总日数
     * @param hd
     * @param rizhao
     * @param h
     * @param lat
     * @return
     */
    public static double caleEtref(double ta, double tmax, double tmin, double wd, double day, double y, double hd, double rizhao, double h, double lat) {
        //湿度换算百分比
        if (hd > 1) {
            hd = hd / 100;
        }
        double pi = 3.141592654;
        double hu = 4.87 / (Math.log(67.8 * 10 - 5.42));
        wd = 4.87 / (Math.log(67.8 * 10 - 5.42)) * wd;
        double p = 101.3 * Math.pow(((273 + ta - 0.0065 * h) / (273 + ta)), 5.26);
        double lam = 2.501 - 0.002361 * ta;

        double r = 0.00163 * p / lam;
        double es = 0.611 * Math.exp(17.27 * ta / (237.3 + ta));
        double slop = 4098 * es /Math.pow( (237.3 + ta), 2);
        double gsc = 0.0820;
        double dr = 1 + 0.033 * Math.cos(2 * day * pi / y);
        double solar = 0.409 * Math.sin(2 * day * pi / y - 1.39);
        lat = pi / 180 * lat;
        double x = Math.max(0.00001, 1 - Math.pow(Math.tan(lat), 2) * Math.pow(Math.tan(solar), 2));
        double ws = pi / 2 - Math.atan(-Math.tan(lat) * Math.tan(solar) / Math.sqrt(x));
        double ra = 24 * 60 / pi * gsc * dr * (ws * Math.sin(lat) * Math.sin(solar) + Math.cos(lat) * Math.cos(solar) * Math.sin(ws));
        double N = 24 / pi * ws;
        double rs = (0.172 + 0.521 * rizhao / N) * ra;
        double rns = 0.77 * rs;
        double rnl1 = 4.903 * Math.pow(10, -9) * (Math.pow(tmax + 273, 4) + Math.pow(tmin + 273, 4)) / 2 * (0.56 - 0.25 * Math.sqrt(es * hd));
        double rnl_xishu = 0.1 + 0.9 * rizhao / N;
        double rnl = rnl1 * rnl_xishu;
        double rn = rns - rnl;

        double etref = (0.408 * slop * rn + r * 900 / (ta + 273) * wd * (es - es * hd)) / (slop + r * (1 + 0.34 * wd));     //!潜在蒸散
//        double etair=(r*900/(ta+273)*wd*(es-es*hd))/(slop+r*(1+0.34*wd));//   !动力项
//        double etrad=(0.408*slop*rn)/(slop+r*(1+0.34*wd));// !辐射项
        return etref;
    }

}
